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Abstract 

In this paper, the process of the introduction of numerical 
diffusion is investigated when low-order differential scheme 
in convection-diffusion equation for discrete suspended load 
is in use. According to the results of instances calculation, 
the optimal choice for parameters when physical diffusion is 
simulated by numerical diffusion is obtained. Finally, 
comparison between several commonly used low-order 
schemes is conducted and the most viable low-order scheme 
in simulating physical diffusion by numerical diffusion is 
proposed in dealing with suspended load. 
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Introduction 

In general, using the low-order scheme for numerical 
calculations will more or less introduce numerical 
diffusion (Xie 1990). Different schemes have varying 
degrees of numerical diffusion. In theory this is not 
conducive to the exact solution of numerical 
calculation, but utilizing the numerical diffusion, such 
as to simulate the physical diffusion phenomenon, is 
also an effective method (Yang 1993). As for how to 
use the numerical diffusion to simulate the physical 
diffusion, and how to select parameters to achieve an 
acceptable level of the simulation results, they are the 
issues to be discussed in this paper. 

Convection-Diffusion Equation for One- 
Dimensional Suspended Load 

Equation for suspended load movement is (Zheng and 
Zhao 2001): 



Introducing the flow continuity equation, the equation 
above can be rewritten as: 

dS QdS 1 d .„ dS (S-r)dA, 

— +— AD x — = - — s - 

dt A dx Adx x dx A dt 

where D x the longitudinal diffusion coefficient, for 
sediment diffusion, it can be approximated as 

D x = 0.25u*h 

Generally, for transportation issue of suspended load, 
the longitudinal diffusion of suspended load is much 
smaller than the longitudinal conviction of sediment, 
which is often overlooked. If the carrier velocity and 
diffusion coefficient D x are constants without 

considering the source item, the convection-diffusion 
equation will be 

as as „ e 2 s 

+ U = D x 2 

dt dx ' dx 

For this equation, the main numerical difficulty is to 
calculate the convection term because it strictly 
demands the conservation of matter, while the 
differential method for numerical solution often 
cannot achieve this. If the solution to convection- 
diffusion is considered under the premise of a good 
solution to convection-diffusion term, the probability 
of successful numerical solution would be higher. If 
ignoring the diffusion term, we can get the pure 
convection equation 

dS dS . 

— + u — = 
dt dx 
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Establishment of Differential Scheme 

Select Upwind scheme discrete pure convection 
equation, when u > , the differential equation is 
(Courant 1928) 

S^ +1 =s n . -C r (S" -s n H ) 

where C r is Crout number and C r = CAt I Ax 
The differential equation can be reduced to 



-S" 



At 



-+C- 



'.H 



Ax 



= 



The first item in the above equation is the forward 
difference quotient of 8S 1 dt , and the second term is 

backward difference quotient of dS 1 dt , then C = u 
can be drawn. We will exam the compatibility, 
stability and convergence of the scheme as follows. 
Expand and S"_ x on point (j,n) 

S^ = s n + {—) n ,At + -{^4-) n At 2 +... 

1 dt ' 2 dt 2 1 

PC 1 P 2 C 

Take them into the differential equation, and use the 
pure convection equation to get 

as as Ax 2 c(\-c,)d 2 s 

- + u — = - 



dt dx 2At dx 1 

Assume p = ^ r (l-C r ) 
2At 



rp. dS dS n 
Then — + u — = D t 

dt dx 



d 2 S 

dx 2 



Comparing to the pure convection equation, only the 
right term of the equation tends to be zero when the 
space and time step are small enough, which means 
they are compatible. Taking At = C r Ax I C into D n , we 

can get D = ^C{\-C r ) 
2 

where C = u ,0 <C r < 1 , when Ax -> , At -> , and 
D n = . So the difference equations and differential 
equations are compatible, as the stability condition for 
Upwind is nAf/Ax<l . Using Lax equivalence 

theorem again we can get that the solution of 
differential equation converges to the solution of 
differential equations. 



From the above proof it is clear that, Upwind scheme 
disperses a pure convection equation and an equation 
similar to convective diffusion can be obtained. 

2 2 

Although the differentiated items D n d S I dx can 

make the equation compatible, calculation shows that 
the item also makes the solution of differential 
equations no longer converge to the solution of the 
original differential equation. Instead it converges to 
the solution of the convection-diffusion equation, 
which means the implicit numerical diffusion, 
resulting from the limitation of Ax , At , which cannot 
be infinitely close to zero. However, since it is known 
that difference equations converge to the convection- 
diffusion equation, then it is used to identify the 
required parameter values, so that the differential 
equations can converge precisely to differential 
equations, namely numerical diffusion can simulate 
the physical diffusion well. 

The Instance Analysis of Parameters Effect 

As differential equation is similar in form with 
convection-diffusion equation, D n might be called 
numerical diffusion coefficient. It is clear that the 
closeness between numerical diffusion coefficient D n 

and physical diffusion coefficient D x can determine 

the accuracy of the numerical diffusion in simulating 
physical diffusion, thus the problem is to find the right 
parameters, so that D n can be close to D ( (Papadakis 
and Metaxas 2011). On the other hand, three 
parameters can determine D which are velocity u , 

space interval Ax and time step At , where flow rate 
should be in accordance with the actual engineering 
value, and Ax i s subject to the restrictions on the 
boundary simulation accuracy (it will take on different 
values in accordance with the different size of the 
project and the complexity of the border. Usually, for 
example, the smaller the size of Ax is, the more 
precise the boundary simulation results will be. But it 
will increase the amount of calculation), only the time 
step At has wider choice. By selecting different At 
we can get the corresponding numerical diffusion 
coefficient D n . Then comparing D n with different 

longitudinal diffusion coefficient D x , the impact of 

At on calculated results can be examined 
(Gasiorowski 2013). 

Build the following model: 
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Model 1 

Length is 10km, width is 2m, water depth is 4m, 
Chezy coefficient C is 25, and flow rate is 0.5m/s, 
Ax =200m. Take respectively At =400s, 399.5s, 399s, 
360s, 300s, 200s... 



The calculation results are shown as follows (Table 1). 

table 1 Calculation results (u=0.5m/s) 



At 


Parameters 


(s) 


c r 


D n 


° x 


400 


1 





0.063246 


399.5 


0.99875 


0.0625 


0.063246 


399 


0.9975 


0.125 


0.063246 


360 


0.9 


5 


0.063246 


300 


0.75 


12.5 


0.063246 


200 


0.5 


25 


0.063246 


100 


0.25 


37.5 


0.063246 


50 


0.125 


43.75 


0.063246 


0.5 


0.00125 


49.9375 


0.063246 



In the above table, only when At =399.5s, the 
difference between numerical diffusion coefficient D n 

and longitudinal diffusion coefficient D f is the 
smallest at 0.000746. At the same time, the difference 
between C r and 1 is 0.00125, when At decreases by 
0.5s to 399s. The two diffusion coefficients differs in an 
order of magnitude, indicating that only when C r is 

very close to 1 the physical diffusion can be simulated 
more accurately. On the other hand, the numerical 
diffusion coefficient D n will increase with the 

decrease of C r . When C r tends 0, a solution with 
larger error will be got. 

Here the time step is accurate to the extent of 0.5, 
which is not conducive to calculation. Since the time 
step may be influenced by the flow velocity, the flow 
rate is changed and the model is modified as follows: 

Model 2 

Flow velocity u =lm/s, taking At =400s, 300 s 200s, 
199.5s, 199s, 198s... other parameters are unchanged. 



In the results (Table 2) when At is greater than Ax, 
numerical diffusion coefficient is negative. It will 
repeat the last calculation law when they are equal. In 
fact if the formula is changed for D n , we can find the 
contradictions. 

_ Ax 2 C r (\-C r ) _ uAx-uAt 
2At ~ 2 

When u =lm/s, D n = (Ax - At) 1 2 , obviously Ax needs 
to be larger than At to make D n positive, in fact, this 
is caused by the failure to satisfy the stability condition 
|k|a*/Ajc<1 . And, from another perspective, the 

accuracy of D n shows direct linear relationship with 

the precision of Ax as well as At , which demands 
highly for Ax and At , in the actual calculation, it is 
often impossible to meet such a request, so other flow 
rates should be taken to calculate again. 



table 2 Calculation results (u=1m/s) 



At 


Parameters 


(s) 


c r 


D n 




400 


2 


-100 


0.126491 


300 


1.5 


-50 


0.126491 


200 


1 





0.126491 


199.5 


0.9975 


0.25 


0.126491 


199 


0.995 


0.5 


0.126491 


198 


0.99 


1 


0.126491 


100 


0.5 


50 


0.126491 


50 


0.25 


75 


0.126491 


0.5 


0.0025 


99.75 


0.126491 



Model 3 

Flow velocity u =0.8m/s, taking Af =250s, 249.5s, 249s, 
200s, 100s, 50s... The remaining parameters are 
unchanged. 

Model 4 

Flow rate u =0.2m/s, taking At =1000s, 999.5s, 999s, 
998s, 900s, 800s... The remaining parameters are 
unchanged. 

Calculation results are shown in Table 3 and Table 4 
respectively. 
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table 3 Calculation results (u=0.8m/s) 



At 


Parameters 


(s) 




D n 


D x 


250 


1 





0.101193 


249.5 


0.998 


0.16 


0.101193 


249 


0.996 


0.32 


0.101193 


200 


0.8 


16 


0.101193 


100 


0.4 


48 


0.101193 


50 


0.2 


64 


0.101193 


0.5 


0.002 


79.84 


0.101193 



table 4 Calculation results (u=0.2m/s) 



At 


Parameters 


(s) 




D n 


D x 


1000 


1 





0.025298 


999.5 


0.9995 


0.01 


0.025298 


999 


0.999 


0.02 


0.025298 


998 


0.998 


0.04 


0.025298 


900 


0.9 


2 


0.025298 


800 


0.8 


4 


0.025298 


500 


0.5 


10 


0.025298 


100 


0.1 


18 


0.025298 


0.5 


0.0005 


19.99 


0.025298 



Combining the above four models; we can draw the 
following rules: 

■ When u is in the interval (0,1]/ the longitudinal 
diffusion coefficient D t increases as u 

increase, and the maximum value is close to 
0.126 (this is the maximum value when the 
water depth is 4m and Chezy coefficient is 25, 
depending on the circumstances, the value will 
be different); 

■ When u is in the interval (0,1], if u is small, 
the numerical diffusion simulation of physical 
diffusion is better, which is the closest when 



u =0.5m/s; 

When u is in the interval (0,1], the simulation 
is better when the time step value is 0.5s 
smaller than the " At that makes C r =1", but 
the difference of 0.5s is not ideal; 

When u is in the interval (0,1], although the 
closer C r is to 1, the better the simulation is, 
however within a small range close to 1, it is 
not in line with this rule. As shown in the 
following table. 

table 5 Variations in parameters with different u 



u 


At 


Parameters 


(m/s) 


(s) 




D n 


D x 


D-D 

x n 


% 


0.2 


999.5 


0.9995 


0.01 


0.025298 


0.0153 


60.48 


0.5 


399.5 


0.99875 


0.0625 


0.063246 


0.0007 


1.11 


0.8 


249.5 


0.998 


0.16 


0.101193 


-0.0588 


58.11 


1.0 


199.5 


0.9975 


0.25 


0.126491 


-0.1235 


97.64 



2.5 - 








O u=0. 2m/s 


2 - 




■ u=0. 5m/ s 






A u=0. 8m/s 


1.5 - 




Xu=l. Om/s 






FIG. 1 Distribution of C r with growth of the time step At (in log 
form) 




FIG. 2 Distribution of D n with growth of the time step At (in log 
form) 

Obviously, in the descending process of u , C r is 
getting closer to 1, but D - D is then minimum when 
u =0.5m/s. The existence of a smaller difference can be 
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tested by changing the value of u , however it is no 
longer studied here. 

In fact, in order to get the parameters conditions of 
numerical diffusion that accurately simulate physical 
diffusion, formulas are used to reason reversely, 
assume D x = D n 
We get: 

^ uAx-u 2 At ^ „ „_ , I — u . 

D n = = D = Q.25u,h = 0.25 Jg —h 

2 v C 

Simplify it to be: 

, O.sJgh 

Ax = uAt + 

C 

By assuming £ _ , we can obtain: 

C 

Ax-k 

At = 

u 

Obviously, At is a linear function of Ax , where k is a 
constant related to the parameters of the physical 
model. Take U =0.5, h =4, C =25, Ax =200m into the 
formula then the calculation shows At =399.499 which 
is very close to the above 399.5, so the difference 
between the two diffusion coefficients is insignificant. 
But this step size is very unfavorable to calculation; as 
it is tendency to take an integer as much as possible. 
Although in different engineering cases, different 
parameters of the physical model can improve the 
value of At , its value depends completely on the 
differential scheme itself. In order to keep the 
generality, the following will horizontally contrast a 
variety of commonly used differential schemes under 



the current model parameter conditions to examine 
whether a scheme can achieve a relatively satisfying 
step value. 

Comparison among Various Difference 
Schemes 

Although various differential schemes have different 
compatibility and stability conditions, as long as 
extracting numerical diffusion coefficient after 
dispersing pure convection equations respectively, 
then equivalent unify with the physical diffusion 
coefficient, we can obtain the calculation formula for 
step size (Table 6). 

Obviously, these three schemes of Dobbins, Upwind 
and implicit have the same method to select the step 
size, and they also have the same characteristics 
represented by the Upwind scheme. The low-order 
Preissmann scheme is more suitable for the large 
stride length, and the possibility of an integral step 
length has been improved in this way, where $7=0.6, 

6 =0.4, At =19900s, and D n =0.05, which shows a 

difference of 20.94% with D x . When At =19870s and 

D n =0.065, D n show a difference of 2.77% with D x 

and precision can meet the requirements. Instead, Lac- 
Friedrichs scheme is more suitable for small time step. 
Though the step is smaller, the simulation accuracy is 
not improved, such as when At =247s and D n =0.097, 

At is only 0.1s smaller than 247.1s, while the relative 
error of diffusion coefficient reaches an unsatisfactory 
level, 52.83%. 



table 6 Results of different calculation formulations 



Scheme name 


D expression 


At expression 


Calculation results 


Dobbins 


Ax 2 C r (\-C r ) 
2At 


Ax-k 

At = 

u 


399.5 


Upwind 


Ax 2 C r {\-C r ) 
2At 


Ax-k 

At = 

u 


399.5 


Implicit 


Ax 2 C r (\ + C r ) 
2At 


k - Ax 

At = 

u 


-399.5 


Preissmann 


2At 


Af feAx + (l-2^)Ax 2 
2(20-1) 


19874.8® 
20125.2® 


Lac-Friedrichs 


Ax 2 (\-C r -C 2 r ) 
2At 


V5Ax 2 +2kAx-k 2 -k-Ax 

At = 

2u 


247.1 



Note: 

1. In the table k = O.S^j(gh) I C , where h =4m, Chezy coefficient C =25, flow rate U =0.5m/s, Ax =200m; 

2. Preissmann scheme only consider the low-order condition, i.e. <p , 9 is not equal to 0.5, ®, ®, are respectively corresponding to the 
different values of tp , 6 , where © is corresponding to p =0.6, and 9 =0.4, © is corresponding to (p =0.4, and 9 =0.6. 
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Conclusions 

The paper starts from Upwind scheme, stating the 
process that generates numerical diffusion when 
dispersing suspended load convection-diffusion 
equation with a low-level difference scheme, and it 
tries to simulate the physical diffusion phenomenon 
using numerical diffusion, simulation accuracy under 
different parameters, such as At , Ax and u , are 
calculated. And by horizontally comparing different 
differential schemes we can get the following 
conclusions: 

1) Numerical diffusion of differential scheme can be 
used to simulate the actual physical diffusion 
phenomenon, but it has certain requirements on 
the parameter selection. For suspended load 
convection-diffusion equation, the simulation 
method is to disperse the pure convection 
equation first, and then propose diffusion term. 
Since the coefficient in front of diffusion term is 
the numerical diffusion coefficient, by assuming 
which equals to the physical diffusion coefficient, 
we can get the relationship equation composed of 
some parameters that meet the simulation 
accuracy. 

2) In general, in the various parameters, only spatial 
interval Ax and time step At which are related 
to the calculation can be chosen randomly, and 
the selection of Ax needs to be compatible with 
the size of the project and the complexity of the 
boundary, so the key is the time step At . By 
calculation, this paper has found the relationship 
of At that meets the accurate simulation under 
various differential schemes, but usually the 
calculated step sizes are with decimals, which is 
not conducive to calculation, so an instance is 
used to calculate the relative error between the 
numerical diffusion coefficient and the physical 
diffusion coefficient when At takes a similar 
integer value. It is found that as the step of low- 
level Preissmann scheme is larger, the possibility 
of an integral step is greatly increased, thereby 



increasing the calculation accuracy, which is a 
recommended method. 

3) Though this article only analyzed suspended load 
convection diffusion equation, it can be 
speculated that bedload convection-diffusion 
equation, as well as the convection-diffusion 
phenomena of pollutants in the water, can be 
analyzed using the same method. Similarly, the 
method can be employed to solve the general 
problems by means of its unique superiority of 
differential equations - clear mathematical 
foundation, simple calculation and easy 
programming, which is still favourably compared 
to a more precise numerical method. 
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